{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "
\n", "
\n", "
Computational Seismology
\n", "
Finite Differences Method with Optimal Operator - Acoustic 1D Case
\n", "
\n", "
\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "

\n", "\n", "\n", "\n", "

\n", "\n", "\n", "---\n", "\n", "This notebook is part of the supplementary material \n", "to [Computational Seismology: A Practical Introduction](https://global.oup.com/academic/product/computational-seismology-9780198717416?cc=de&lang=en&#), \n", "Oxford University Press, 2016.\n", "\n", "##### Authors:\n", "* Heiner Igel ([@heinerigel](https://github.com/heinerigel))\n", "* Lion Krischer ([@krischer](https://github.com/krischer))\n", "* Taufiqurrahman ([@git-taufiqurrahman](https://github.com/git-taufiqurrahman))\n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook covers the following aspects:\n", "\n", "* implementation of the 1D acoustic wave equation\n", "* understanding the input parameters for the simulation and the plots that are generated\n", "* allowing you to explore the finite-difference method with optimal operator\n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Optimising Operators\n", "\n", "Geller and Takeuchi (1995) developed criteria against which the accuracy of frequency-domain calculation of synthetic seismograms could be optimised. This approach was transferred to the time-domain finite-difference method for homogeneous and heterogeneous schemes by Geller and Takeuchi (1998). Look at an optimal operator and compare with the classic scheme. The space-time stencils are illustrated below\n", "\n", "\n", "\n", "$$Conventional-(1/dt^2)$$\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
t+dt1
t-2
t-dt1
x-dxxx+dx
\n", "\n", "$$Optimal-(1/dt^2)$$\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
t+dt1/1210/121/12
t-2/12-20/12-2/12
t-dt1/1210/121/12
x-dxxx+dx
\n", "\n", "\n", "$$Conventional-(1/dx^2)$$\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
t+dt
t1-21
t-dt
x-dxxx+dx
\n", "\n", "$$Optimal-(1/dx^2)$$\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
t+dt1/12-2/121/12
t10/12-20/1210/12
t-dt1/12-2/121/12
x-dxxx+dx
\n", "\n", "*The conventional 2nd order finite-difference operators for the 2nd derivative are compared with the optimal operators developed by Geller and Takeuchi (1998), see text for details.*\n", "\n", "\n", "Note that summing up the optimal operators one obtains the conventional operators. This can be interpreted as a smearing out of the conventional operators in space and time. The optimal operators lead to a locally implicit scheme, as the future of the system at $(x,t+dt)$ depends on values at time level \"$t+dt$, i.e., the future depends on the future. That sounds impossible, but it can be fixed by using a predictor-corrector scheme based on the first-order Born approximation. \n", "\n", "\n", "The optimal operators perform in a quite spectacular way. With very few extra floating point operations an accuracy improvement of almost an order of magnitude can be obtained. The optimal scheme performs substantially better than the conventional scheme with a 5-point operator.\n", "\n", "---" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "# Import Libraries (PLEASE RUN THIS CODE FIRST!) \n", "# ----------------------------------------------\n", "import numpy as np\n", "import matplotlib\n", "# Show Plot in The Notebook\n", "matplotlib.use(\"nbagg\")\n", "import matplotlib.pyplot as plt\n", "\n", "# Sub-plot Configuration\n", "# ----------------------\n", "from matplotlib import gridspec \n", "\n", "# Ignore Warning Messages\n", "# -----------------------\n", "import warnings\n", "warnings.filterwarnings(\"ignore\")" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "# Parameter Configuration \n", "# -----------------------\n", "nt = 501 # number of time steps\n", "eps = 1. # stability limit\n", "xs = 250 # source location in grid in x-direction\n", "xr = 450 # receiver location in grid in x-direction\n", "\n", "# Material Parameters\n", "# -------------------\n", "rho = 2500. # density\n", "c0 = 2000. # velocity\n", "mu = rho*(c0**2) # elastic modulus \n", "\n", "# Space Domain\n", "# ------------\n", "nx = 2 *xs # number of grid points in x-direction\n", "dx = 2.*nx/(nx) # calculate space increment\n", "\n", "# Calculate Time Step from Stability Criterion\n", "# --------------------------------------------\n", "dt = .5*eps*dx/c0" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "# Source Time Function \n", "# --------------------\n", "f0 = 1./(10.*dt) # dominant frequency of the source (Hz)\n", "t0 = 4./f0 # source time shift\n", "\n", "# Source Time Function (Gaussian)\n", "# -------------------------------\n", "src = np.zeros(nt)\n", "time = np.linspace(0 * dt, nt * dt, nt)\n", "# 1st derivative of a Gaussian\n", "src = -2.*(time-t0)*(f0**2)*(np.exp(-1.*(f0**2)*(time-t0)**2))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "# Operator Error Calculation \n", "# --------------------------\n", "\n", "# Conventional FD Operators\n", "# -------------------------\n", "A0 = rho / (dt**2) * np.matrix\\\n", "(' 0. 1. 0.;\\\n", " 0. -2. 0.;\\\n", " 0. 1. 0.')\n", "K0 = mu / (dx**2) * np.matrix\\\n", "(' 0. 0. 0.;\\\n", " 1. -2. 1.;\\\n", " 0. 0. 0.')\n", "\n", "# Modified FD Operators\n", "# ---------------------\n", "A = 1. / 12. * rho / (dt**2) * np.matrix\\\n", "(' 1. 10. 1.;\\\n", " -2. -20. -2.;\\\n", " 1. 10. 1.')\n", "K = 1. / 12. * mu / (dx**2) * np.matrix\\\n", "(' 1. -2. 1.;\\\n", " 10. -20. 10.;\\\n", " 1. -2. 1.')\n", "\n", "# Calculate Operator Error\n", "# ------------------------\n", "dA = A0 - A # error of conventional operator A\n", "dK = K0 - K # error of conventional operator K\n", "d0 = (dA - dK) * (dt**2) / rho # basic error of modified operator" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "code_folding": [] }, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support. ' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
');\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('